#include <cusp/array2d.h>
#include <cusp/multiply.h>
#include <cusp/print.h>
#include <cusp/qr.h>

#include <limits>
#include <iostream>


int main(void)
{
	typedef typename cusp::host_memory MemorySpace;
	typedef float ValueType;

	cusp::array2d<ValueType, MemorySpace> A(3,3);
	A(0,0) = 12.0; A(0,1) = -51.0;  A(0,2) = 4.0;
	A(1,0) =  6.0; A(1,1) = 167.0;  A(1,2) = -68.0;
	A(2,0) = -4.0; A(2,1) =  24.0;  A(2,2) = -41.0;
	
	cusp::array2d<ValueType, MemorySpace> Q, R;

	std::cout << "Matrix A : ";
	cusp::print(A);
	std::cout << "Matrix R : ";
	cusp::print(R);

	std::cout << "Perfoming QR decomposition..." << std::endl;
	cusp::qr(A, Q, R);

	std::cout << std::setprecision(std::numeric_limits<ValueType>::digits10 + 1);

	std::cout << "Matrix A : ";
	cusp::print(A);
	std::cout << "Matrix Q : ";
	cusp::print(Q);
	std::cout << "Matrix R : ";
	cusp::print(R);

	cusp::multiply(Q, R, A);
	std::cout << "Matrix Q*R : ";
	cusp::print(A);
	

	return 0;
}

